rm(list=ls())
############ initialize  ##########################
setwd("/Volumes/Disk2/Future_PM/Dataverse")
library(fields); library(maps); library(ncdf4);library(abind)
source('Function/get_geo.R')
source('Function/get_met.R')
source('Function/functions.R')

zz<-read.table("Data/US_Mask.txt",header=T)
zz=t(zz)
zz=zz[,(dim(zz)[2]):1]
region.ind=(zz>=1)
US_region=zz>0

load("Data/PM_region.Rdata")
US_region=PM.region==1

wd=getwd()
########################### SEASON #############################
cal.season=function(start_month,end_month){
filenames= Sys.glob(paste(wd,'/Data/Future_PM25/*.Rdata',sep=""))
for (dir.num in 1:length(filenames)){	
	print(dir.num)
	ss=load(filenames[dir.num])
	CMIP5.hybrid.PM[is.infinite(CMIP5.hybrid.PM)]=NA
	if(start_month<=end_month){
	    ind1=(CMIP5.time[,1]<=2019 & (CMIP5.time[,2]>= start_month & CMIP5.time[,2]<= end_month))
	    ind2=(CMIP5.time[,1]>=2050 & (CMIP5.time[,2]>= start_month & CMIP5.time[,2]<= end_month))
    }else {
	    ind1=(CMIP5.time[,1]<=2019 & (CMIP5.time[,2]>= start_month | CMIP5.time[,2]<= end_month))
	    ind2=(CMIP5.time[,1]>=2050 & (CMIP5.time[,2]>= start_month | CMIP5.time[,2]<= end_month))    	
    }
     change1=apply(CMIP5.hybrid.PM[,,ind2],c(1,2),mean,na.rm=TRUE)-apply(CMIP5.hybrid.PM[,,ind1],c(1,2),mean,na.rm=TRUE)	
 	if (dir.num==1) {
		total_change1=change1
	}
	if (dir.num>1)  {
		total_change1=abind(total_change1,change1,along=3)
	}
}
return(total_change1)
}

rwb.colors=tim.colors
######### figures
dev.new(width=6,height=4)
mai=c(0.1, 0.05, 0.2, 0.05)
mgp=c(1.2, 0.4, 0)
tcl=-0.2
ps=12


m <- rbind(c(0,1, 0.93, 1),
c(0, 0.33, 0.6, 0.95),c(0.33, 0.67,0.6,0.95),c(0.67, 1,0.6,0.95), 
c(0, 0.33, 0.25, 0.6),c(0.33, 0.67,0.25,0.6),c(0.67, 1,0.25+0.1,0.6+0.1))
close.screen(all.screens = TRUE)
split.screen(m)

screen(1)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
text(x=0.5,y=0.5,bquote('Effects of 2050s climate change on '*PM[2.5]~ ''),cex=1.1)

screen(2)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
spdata=cal.season(3,5)
ap=apply(spdata,c(1,2),mean,na.rm=TRUE)
ap[! US_region]=NA
num=apply(spdata>0,c(1,2),sum,na.rm=TRUE);sig=(num>=14 | num <= 5)
ap[!sig]=NA
image(lon.frame,lat.frame,ap,zlim=c(-3.5,3.5),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='')
ind=which(sig & PM.region==1,arr.ind=TRUE)
# points(lon.frame[ind[,1]],lat.frame[ind[,2]],pch=16,cex=0.2)
map("world",add=TRUE)
title('(a) MAM',font.main=1,cex.main=1)

screen(3)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
spdata=cal.season(6,8)
ap=apply(spdata,c(1,2),mean,na.rm=TRUE)
ap[! US_region]=NA
num=apply(spdata>0,c(1,2),sum,na.rm=TRUE);sig=(num>=14 | num <= 5)
ap[!sig]=NA
image(lon.frame,lat.frame,ap,zlim=c(-3.5,3.5),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='')
ind=which(sig & PM.region==1,arr.ind=TRUE)
# points(lon.frame[ind[,1]],lat.frame[ind[,2]],pch=16,cex=0.2)
map("world",add=TRUE)
title('(b) JJA',font.main=1,cex.main=1)

screen(4)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
spdata=cal.season(9,11)
ap=apply(spdata,c(1,2),mean,na.rm=TRUE)
ap[! US_region]=NA
num=apply(spdata>0,c(1,2),sum,na.rm=TRUE);sig=(num>=14 | num <= 5)
ap[!sig]=NA
image(lon.frame,lat.frame,ap,zlim=c(-3.5,3.5),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='')
ind=which(sig & PM.region==1,arr.ind=TRUE)
# points(lon.frame[ind[,1]],lat.frame[ind[,2]],pch=16,cex=0.2)
map("world",add=TRUE)
title('(c) SON',font.main=1,cex.main=1)

screen(5)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
spdata=cal.season(11,2)
ap=apply(spdata,c(1,2),mean,na.rm=TRUE)
ap[! US_region]=NA
num=apply(spdata>0,c(1,2),sum,na.rm=TRUE);sig=(num>=14 | num <= 5)
ap[!sig]=NA
image(lon.frame,lat.frame,ap,zlim=c(-3.5,3.5),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='')
ind=which(sig & PM.region==1,arr.ind=TRUE)
map("world",add=TRUE)
title('(d) DJF',font.main=1,cex.main=1)

screen(6)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
spdata=cal.season(1,12)
ap=apply(spdata,c(1,2),mean,na.rm=TRUE)
ap[! US_region]=NA
num=apply(spdata>0,c(1,2),sum,na.rm=TRUE);sig=(num>=14 | num <= 5)
ap[!sig]=NA
image(lon.frame,lat.frame,ap,zlim=c(-3.5,3.5),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='')
ind=which(sig & PM.region==1,arr.ind=TRUE)
map("world",add=TRUE)
title('(e) Annual',font.main=1,cex.main=1)

screen(7)
 mai=c(0.1, 0.1, 0.05, 0.2)
par(mai=mai, mgp=mgp, tcl=tcl, ps=8.5)
image.plot(legend.shrink=1,legend.only=TRUE,zlim=c(-3.5,3.5),col=rwb.colors(32),horizontal=TRUE,legend.width=1.8,axis.args = list(mgp = c(0, 0.5, 0),tcl=-0.2,padj=-1
))
text(bquote('unit: '~mu*'g' ~m^-3*''),x=par("usr")[2]-0.11,y=par("usr")[3]+0.31,srt=0, adj = 0,xpd=TRUE,cex=1.3)

